# define objects to save model estimates for plotting
diffs <- NA
diffs.se <- NA
beta1 <- beta2 <- NA
se1 <- se2 <- NA


VDEM <- read.csv("V-Dem-DS-CY-v6.1.csv")
VDEM <- subset(VDEM, select=c(COWcode, year, v2cltort, v2cltort_sd, v2clkill, v2clkill_sd, country_name))
names(VDEM) <- c("COW", "YEAR", "v2cltort", "v2cltort_sd", "v2clkill", "v2clkill_sd", "country_name")

year <- 1949:2010


# loop through all the different model specifications for this treaty variable
for(j in 1:length(year)){

data <- VDEM

OBS <- rev(cumsum(rev(as.numeric(table(data$YEAR)))))

# load binary treaty variables
treaty <- read.csv("hr_treaties.csv")
treaty <- subset(treaty, select=c(ccode, year, cat_rn, cescr_rn, cerd_rn, ccpr_rn, cedaw_rn, crc_rn))
names(treaty) <- c("COW", "YEAR", "cat", "cescr", "cerd", "ccpr", "cedaw", "crc")

# this is not correct but necessary for the replications
WRONG <- TRUE
if(WRONG==TRUE)treaty$cat[is.na(treaty$cat)] <- 0
if(WRONG==TRUE)treaty$cescr[is.na(treaty$cescr)] <- 0
if(WRONG==TRUE)treaty$cerd[is.na(treaty$cerd)] <- 0
if(WRONG==TRUE)treaty$ccpr[is.na(treaty$ccpr)] <- 0
if(WRONG==TRUE)treaty$cedaw[is.na(treaty$cedaw)] <- 0
if(WRONG==TRUE)treaty$crc[is.na(treaty$crc)] <- 0

# load polity data
polity <- read.csv("p4v2010.csv", na.strings = c("-66", "-77", "-88", "-99"))
polity2 <- subset(polity, select=c(ccode, year, polity2))

# load updated Gleditsch trade gdp and population data
new.gdp <- read.delim("gdpv6.txt")
new.gdp <- subset(new.gdp, select=c(statenum, year, rgdppc, pop))
names(new.gdp) <- c("ccode", "year", "rgdpch", "pop")

# merge the data together and make lags
data <- merge(data, treaty, by.x=c("COW", "YEAR"), by.y=c("COW", "YEAR"), all.x=TRUE, all.y=FALSE, suffixes = c("",".treaty"))
nrow(data)

data <- merge(data, polity2, by.x=c("COW", "YEAR"), by.y=c("ccode", "year"), all.x=TRUE, all.y=FALSE)
nrow(data)

data <- merge(data, new.gdp, by.x=c("COW", "YEAR"), by.y=c("ccode", "year"), all.x=TRUE, all.y=FALSE)
nrow(data)


data.lag <- data
data.lag$YEAR <- data.lag$YEAR + 1
data <- merge(data, data.lag, by.x=c("COW", "YEAR"), by.y=c("COW", "YEAR"), all.x=TRUE, all.y=FALSE, suffixes = c("",".lag"))

data <- subset(data, !is.na(v2cltort.lag) & !is.na(cat.lag) & !is.na(polity2.lag) & !is.na(rgdpch.lag) & !is.na(pop.lag) & YEAR>=year[j] & YEAR<=2010)
nrow(data)

newdata <- list()

# take n draws from the posterior distribution and make n new datasets in the list object just created
for(i in 1:50){
  newdata[[i]] <- data
  newdata[[i]]$draw.lag <- rnorm(cbind(rep(1,nrow(data))), mean=data$v2cltort.lag, sd=data$v2cltort_sd.lag)
  newdata[[i]]$draw <- rnorm(cbind(rep(1,nrow(data))), mean=data$v2cltort, sd=data$v2cltort_sd)
}


# define model formula for each specification
#FML <- "v2cltort ~ draw.lag + cat.lag"
#FML <- "v2cltort ~ draw.lag + cat.lag + polity2.lag"
#FML <- "v2cltort ~ draw.lag + cat.lag + polity2.lag + log(rgdpch.lag)"
FML <- "v2cltort ~ draw.lag + cat.lag + polity2.lag + log(rgdpch.lag) + log(pop.lag)"
#FML <- "v2cltort ~ draw.lag + cat.lag + log(rgdpch.lag) + log(pop.lag)"
#FML <- "v2cltort ~ draw.lag + cat.lag + log(rgdpch.lag)"
#FML <- "v2cltort ~ draw.lag + cat.lag + log(pop.lag)"
#FML <- "v2cltort ~ draw.lag + cat.lag + polity2.lag + log(pop.lag)"

# define function to combine multe model estimates (this incorporates uncertatiny from the latent variables, see the supplementary appendix section F for more details about this function)
milm <- function(fml, midata){
 xx <- terms(as.formula(fml))
 lms <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
 ses <- matrix(data=NA, nrow=(length(attr(xx, "term.labels")) + 1), ncol=length(midata))
 vcovs <- list()
  for(i in 1:length(midata)){
    tmp <- lm(formula=as.formula(fml), data=midata[[i]])
    lms[,i] <- tmp$coefficients
    ses[,i] <- sqrt(diag(vcov(tmp)))
    vcovs[[i]] <- vcov(tmp)
  }
 par.est <- apply(lms, 1, mean)
 se.within <- apply(ses, 1, mean)
 se.between <- apply(lms, 1, var)
 se.est <- sqrt(se.within^2 + se.between*(1 + (1/length(midata))))
 list("terms"=names(tmp$coefficients), "beta" = par.est, "SE"=se.est, "vcovs"=vcovs,"coefs" = lms)    
}

# call the milm function with the list of datasets newdata
model <- milm(fml=FML, midata=newdata)

# create object for printing during the loop and then print the results to screen
temp <- data.frame(model$terms, model$beta, model$SE, model$beta/model$SE)

print(paste("Model Output"))
print(temp)

print(as.numeric(temp[temp$model.terms=="cat.lag",2]))
print(as.numeric(temp[temp$model.terms=="cat.lag",3]))



#par(mar=c(5,7,5,5))

VDEM.est <- as.numeric(temp[temp$model.terms=="cat.lag",2])
VDEM.se <- as.numeric(temp[temp$model.terms=="cat.lag",3])



beta1[j] <- VDEM.est
se1[j] <- VDEM.se

}





# barplot of the differences across all 8 model specifications
par(mar=c(4,4,4,6))
par(mfrow=c(1,1))
barplot(beta1, space=0, ylim=c(-.15,.25), col=grey(.9), yaxt="n", border="white", xlim =c(1,60))
abline(h=0, col=grey(.5))
for(i in 1:7){abline(v=(c(2,12,22,32,42,52,62)-.5)[i], lwd=.75, col=grey(.75), lty=2)}
for(i in 1:length(year)){
    lines(c(i-.5,i-.5), c(beta1[i]+1.96*se1[i], beta1[i]-1.96*se1[i]), lwd=1)
    lines(c(i-.5,i-.5), c(beta1[i]+se1[i], beta1[i]-se1[i]), lwd=2)
}
#box()
#axis(side=1, at=c(.5:61.5), labels=rep("", 62))
axis(2, at=c(-.20, -.15, -.10, -.05, 0, .05, .10, .15, .20), las=2)
mtext(side=4, "Model Coefficient \nChanging Standard DV", line=3.0, cex=1.0)
mtext(side=3, "Human Rights Latent Treaty Variable Coefficients", line=-1.5, cex=1.25)

axis(side=1, at=c(.5:61.5), labels=rep("", 62))
axis(1, at=c(2,12,22,32,42,52,62)-.5, label=c(1950, 1960, 1970, 1980, 1990, 2000, 2010), tick=F, cex.axis=1.5)
axis(2, at=c(-.35,-.30, -.25, -.20, -.15, -.10, -.05, 0, .05, .10, .15, .20, .25, .30, .35), las=2)

